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By using a full many-body approach, we calculate the excitation energy, the effective mass and the 
density profile of soliton states in a three dimensional Bose gas of hard spheres at zero temperature. 
The many-body wave fnnction nsed to describe the soliton contains a one-body term, derived from 
the solution of the Gross-Pitaevskii eqnation, and a two-body Jastrow term which accounts for 
the repulsive correlations between atoms. We optimize the parameters in the many-body wave 
function via a Variational Monte Carlo procedure, calculating the grand-canonical energy and the 
canonical momentum of the system in a moving reference frame where the soliton is stationary. 
As the density of the gas is increased, signihcant deviations from the mean-field predictions are 
found for the excitation energy and the density profile of both dark and grey solitons. In particular, 
the soliton effective mass m* and the mass mAN of missing particles in the region of the density 
depression are smaller than the result from the Gross-Pitaevskii equation, their ratio being however 
well reproduced by this theory up to large values of the gas parameter. We also calculate the profile 
of the condensate density around the soliton notch finding good agreement with the prediction of 
the local density approximation. 

PACS numbers: 67.85.De, 03.75.Kk, 03.75.Lm, 05.30.Jp 


I. INTRODUCTION 

Solitons are nonlinear collective excitations which ap¬ 
pear in a wide range of physical systems, from classical 
fluids to optical fibers: they are characterized as local¬ 
ized wave forms which travel in a uniform medium at 
a constant velocity without spreading. Ultracold atomic 
gases are well-controlled quantum systems which are par¬ 
ticularly suitable for the investigation of solitons. Dark 
and grey solitons (i.e. localized density depressions in 
a homogeneous background) can indeed be produced in 
repulsive Bose-Einstein condensates (BEC^ by imprint¬ 
ing a phase juim on the atomic cloud [iH^, by density 
engineering p, l6| or by merging two coherent BECs ini¬ 
tially prepared in a double well potential Q- Further¬ 
more, other kinds of solitons (e.g. bright solitons, gap 
solitons and dark-bright solitons) have also been created 
and detected in bosonic quantum gases I, [M3. 

From the theoretical point of view, soliton excitations 
in superfluids are studied mostly within mean-field ap¬ 
proaches based on the description of the system in terms 
of a complex order parameter which evolves in space 
and time according to a nonlinear equation. In the case 
of Bose superfluids the paradigmatic theory is provided 
by the Gross-Pitaevskii (GP) equation where soliton so¬ 
lutions have been widely investigated . A concep¬ 
tually similar, even though technically more involved, 
mean-field approach exists also for Fermi superfluids and 
is based on the time-dependent Bogoliubov-de Gennes 
equations (l3 - [T5l| - We point out that mean-field soli¬ 
tons are stable excitations in one dimension (ID), where 
the phase of the order parameter changes sign at a sin¬ 


gle point. In higher dimensions they can be stabilized 
by a sufficiently strong confinement which reduces the 
soliton nodal surface, but otherwise undergo a snake 
instability towards the formation of vortices or vortex 
rings [l,[il[I3|- 

An important question is to understand how the prop¬ 
erties of soliton excitations change when interparticle in¬ 
teractions and/or correlations increase and the mean- 
field picture in terms of the order parameter eventually 
fails. Using time-dependent Bogoliubov theory 11814^ 
as well as more sophisticated numerical techniques [2ll - 
[13 to treat the many-body dynamics of ID bosons at 
zero temperature, it has been shown that beyond mean- 
field effects tend to deplete the condensate and to fill 
the soliton notch making dark solitons unstable. To our 
best knowledge, however, no theoretical studies have ad¬ 
dressed the role of many-body correlations in determining 
relevant properties of grey solitons in a three-dimensional 
(3D) Bose gas, such as their excitation energy, density 
profile and effective mass. In particular, the problem of 
the value of the effective mass beyond the limits of mean- 
field theory became prominent in the first interpretation 
of the recent experiment [13 , where localized excitations 
obtained by phase imprinting in an ultracold Fermi gas 
of ®Li were observed to move much more slowly than pre¬ 
dicted by mean-field theory. These defects were initially 
ascribed as solitons, but lately a more precise imaging 
technique revealed that the moving defects were solitonic 
vortices, thus explaining the puzzle of their greater iner¬ 
tia [13 • Nevertheless, the dependence of the soliton effec¬ 
tive mass on the strength of interactions as the gas gets 
less dilute remains an open question worth considering. 

In this paper we perform full many-body simulations of 


soliton excitations in a 3D Bose gas at zero temperature. 
We calculate the energetic and structural properties of 
grey and dark solitons, both in the very dilute and in 
the strongly interacting regime, and we determine their 
effective mass from the excitation energy vs. speed de¬ 
pendence. The results are systematically compared with 
the predictions of GP equation and important deviations 
are found as the density is increased. Remarkably, the 
ratio between the soliton effective mass and the 

mass of missing particles in the density dip, which is the 
crucial parameter in the dynamics of solitons as macro¬ 
scopic localized defects, is found to remain close to the 
mean-field value up to very high densities. 

The paper is organized as follows. In Sec. m we in¬ 
troduce the variational approach used in the numerical 
determination of the properties of the soliton: we discuss 
the form of the many-body wave function and the energy 
functional in terms of which the soliton state can be de¬ 
fined. In Sec. uni we present the results obtained within 
this many-body approach and we discuss the appearance 
of corrections beyond the mean-field description of soli- 
tons. Finally, in Sec. IlYl we draw our conclusions. 


II. NUMERICAL METHOD 


The system of N interacting identical bosons of mass 
m is described by the Hamiltonian of the quantum de¬ 
generate hard-sphere (HS) model, that is 


H = - 


2m 


i=l i<j 


( 1 ) 


where is the distance between the i-th and j-th par¬ 
ticle and the central potential V is modeled by the hard- 
sphere interaction. 


V{r) 


00 (r < a) 

0 (r > a) , 


( 2 ) 


in terms of the s-wave scattering length a. This model 
has been widely used to describe quantum many-body 
systems with short-range repulsive interaction, both in 
the weakly and in the strongly interacting regime. In¬ 
deed, it is able not only to capture the essential prop¬ 
erties of dilute systems like ultracold gases with positive 
scattering-length, for which the details of the interatomic 
potential are irrelevant, but also to characterize semi- 
quantitatively the static properties of strongly interact¬ 
ing systems like superfluid "‘He [2§ - l^ . Moreover, the 
choice of the HS model is particularly convenient as it 
allows to parametrize the strength of the interactions by 
just varying the value of the dimensionless gas parameter 
na^, where n = N/V is the density of the gas. 

To describe soliton excitations, we consider states of 
the gas described by the many-body wave function 


N 

ff's(ri,..., r^v; t) = Ms (t>{zi - vt) (3) 

i=l i<j 


where As is a normalization factor. The product of 
two-body Jastrow terms accounts for interparticle cor¬ 
relations: we model them using the positive function 
f(r) = in the interval a < r < Rm, that is 

the exact solution of the scattering problem for a hard- 
sphere potential. For r > Rm, we match /(r) with the 
constant . This function satisfies the bound¬ 

ary condition /(r) = 0 at r < a and the wave vector 
k is determined from the condition f'(r = Rm) = 0 of 
the first derivative f at the matching point Rm- The 
value of Rm is used as a variational parameter and is 
optimized to minimize the excitation energy of the dark 
soliton. The complex one-body term cj) is given by 

(l){z) = \/l — tanh 1 — _|_ jq, ^ (qj 

and describes a perturbation in the density profile prop¬ 
agating along the 2 -axis with velocity v. Here ^ = 
l/'/Sirna is the healing length of the gas at the back¬ 
ground density n, while a and 7 are two variational pa¬ 
rameters. The functional form of (j) is dictated from the 
solution ,/n(j){z — of the GP equation describing dark 
and grey solitons [TJ: in this case, the values of the pa¬ 
rameters are 7 = \/2 and a = vIcq, with cn = ^ , the 

' / 'J V2m^ 

speed of sound within the GP theory. The phase of (f) 
undergoes the finite change AS = 2 arccos(a) as 2 varies 
from -boo to — 00 , whereas the density perturbation is 
localized around the origin of the 2 -axis on a length scale 
fixed by y^/Vl — cA ■ The wave function m is transla- 
tionally invariant in the xy plane and its global phase 
is the sum of single-particle contributions each of which 
changes sign at the Zi — vt = 0 plane. 

Since we assumed that the time t enters the many- 
body wave function only through the linear combination 
Zi — vt, it is convenient to describe the problem in a 
moving reference frame, where the soliton is stationary. 
In analogy with the GP equation we define soliton states 
as stationary points of the functional 

= , (5) 

where E is the grand-canonical energy of the many-body 
system, 


E[^s] = J d'R m'R)iH - /i)^s(R) , (6) 


and Pc is the canonical momentum 

iTl f 

-Pelt's] = finLr,Ly{AS - n) - — / 


dR 


( 7 ) 


N , N 


i=l 


i=l 


The set of particle coordinates R = (fi,...,fAr) refers 
to the moving reference frame = (xi,yi,Zi), where 
Zi = Zi — vt. Furthermore, and Ly indicate the length 












of the system respectively in the x and y directions. In 
Eq. (IH), /r is the chemical potential of the homogeneous 
gas enforcing the boundary condition that the soliton 
state reaches the asymptotic density n when Zj = ±oo for 
all the particles. The first term in Eq. 0 arises instead 
from the boundary condition that the many-body wave 
function has a phase jump NAS across the soliton 
and it is necessary to account for the motion of the back¬ 
ground density in the moving reference frame [ill . [l3 . . 
The functional SI can be derived from the principle of 
minimal action applied to the quantity 


A = 


dt 





z/5 


where 4's is a state of the form ([3]) c omp lying with the 
boundary conditions explained above [32|. 

The simulations for a given value of the gas parameter 
na^ are carried out using the Variational Monte Carlo 
(VMC) technique. As a first step, we calculate the prop¬ 
erties of the uniform ground state by performing VMC 
simulations of Njj particles in a box of volume V = Njj/n 
with periodic boundary conditions. We use the func¬ 
tion \l'j/(R) = Au Y\i<j to describe the normalized 

ground state and we determine the corresponding energy 
Eq and chemical potential /i = The soliton state 

is simulated in the same box of volume V = L^LyEz, us¬ 
ing a number Ns < Njj of particles chosen so as to com¬ 
ply with the boundary condition of the density profile 
reaching the uniform value n far away from the soliton 
plane. The volume of the simulation box is chosen large 
enough to prevent finite-size effects: typical values are 
Lz ^ 20^ and L^; = Ly > 8 ^. Moreover, the value Rm 
of the matching point in the Jastrow factor is kept the 
same for the calculation of the and dtg state. For each 
value of the velocity v of the soliton (given as an input 
parameter), we calculate as a function of the varia¬ 
tional parameters and we look for the stationary points 
to determine the optimal values of a and 7 . In analogy 
with the GP equation, the soliton state corresponds to a 
saddle point oi SI, i.e. a minimum as a function of 7 and 
a maximum as a function of a. 


III. RESULTS 

We start the discussion of the structural properties by 
considering the limit of vanishing velocity, u/cq — t 0, 
where the optimal value of a approaches zero and the 
wave function dts describes a dark soliton with a phase 
jump AS = TT. The density profiles n{z) of these station¬ 
ary excitations are shown in Fig. [T] for different values 
of the gas parameter and compared with the profile ob¬ 
tained within the GP approach. A systematic analysis 
of the effects of many-body interactions on the conden¬ 
sate wave function in time-dependent problems has been 
performed in Ref. [s^. In this work, it has been shown 
that the relevant parameter for the deviations from the 


Figure 1. (Color online) Density profile of a dark soliton for 
different values of the gas parameter (lines are a guide to the 
eyes). The solution of the GP equation is also shown for 
comparison. 
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Figure 2. (Color online) Density profile of the soliton at the 
density na^ = 10“^ for different velocities in the reference 
frame where the soliton is at rest. The mean-field curves are 
also shown. 


solution of GP equation is the square root of the quan¬ 
tum depletion, which scales as In Fig.[Tl we can 

see that the VMG results are in excellent agreement with 
mean-field predictions at the lowest density, na^ = I 0 “®. 
As the gas parameter increases, we notice some devia¬ 
tions from the solution of GP equation although they 
become significative only at na^ > 10“^. At these higher 
densities the width of the soliton, in units of the heal¬ 
ing length, decreases and some oscillations appear in the 
density profile. These oscillations are more pronounced 
at the largest densities and get damped in regions far 
from the soliton. Similar oscillations appear in the den¬ 
sity profile of liquid "‘He around a vortexjs^ and their 
appearance can be interpreted as the tendency of the sys¬ 
tem to organize itself in shells of atoms around the defect. 
We notice that, for all values of the gas parameter, the 
density vanishes on the z = 0 soliton plane in agreement 
with the GP prediction for dark solitons. This is a conse- 











Figure 3. (Color online) Condensate density profile no(z) of a 
dark soliton at na® = 10~^. The results of both GP equation 
and LDA are also shown. Inset: local condensate fraction 
naiyZ)!n{z) compared with the LDA result. 


qtience of the choice (0]) of the one-body term which, as 
a —^ 0, generates for each particle a phase discontinuity 
at the z = 0 plane. 

In Fig. m we show the density profile of grey solitons 
moving with different velocities u in a gas of background 
density na^ = 10“^, which corresponds to a regime of 
relatively strong interactions among the particles. One 
observes clear deviations from the GP profile in the re¬ 
gion of the wings of the soliton, similar to the = 0 case 
(see Fig. [T]). The discrepancies between the two profiles 
are less pronounced in the central region of the soliton 
where the VMC results at all velocities are only slightly 
above the GP predictions. 

Another important quantity accessible in our micro¬ 
scopic approach is the condensate profile no{z) around 
the defect. This can be obtained from the long-distance 
behavior of the one-body density matrix in the xy plane 
at a fixed value of the z coordinate: 

2 ^s(r'i,r2,.. ■ ,rjv) 
^'s(ri,r2 ,...,rAf) ’ 

(9) 

with r'l = (a;i,?/i,z), ri = (xi,?/i,z) and R^-i = 
{r2 ,..., Tat}. In Fig. [3] we show the condensate profile 
of a dark soliton (v = 0) at the density na^ = 10“^. Far 
away from the soliton plane at z = 0, the condensate den¬ 
sity reaches the bulk value Uq ~ 0.8n in agreement with 
the result for a homogeneous gas [s^. When approaching 
the soliton plane, we find that no(z) is in closer agree¬ 
ment with the GP profile than the total density n(z). 
Remarkably, the condensate profile is well reproduced by 
the local density approximation (LDA) where we deter¬ 
mine no(z) from the results of the quantum depletion in a 
homogeneous gas [lElal at the local density n(z). Only 
in the region very close to the soliton plane, the LDA re¬ 
sult shows deviations on the order of 5% with respect to 
the local condensate fraction no(z)/n(z) obtained from 
VMG calculations (see inset of Fig. 131). 

The energetics of solitons is reported in Fig. 01 where 
we show the excitation energy AE of the soliton as a 


no(z) = lim 

k'i-ri|->-oo , 


dR^-i |d)'s| 


Figure 4. (Color online) Excitation energy AE of the soliton 
as a function of (v/co)^. The solid line corresponds to the GP 
result. Statistical errors are below symbol size. 
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Figure 5. (Color online) Ratio of the VMC to the GP effective 
mass m* and mass of missing particles in the soliton dip mAN 
as a function of the gas parameter. Inset: ratio as a 

function of na^. 



function of v'^ and for different values of the gas pa¬ 
rameter. This value is obtained from the difference 
AE = Sill's] — E[^!u] of the grand-canonical energy dHl) 
between the states with and without the soliton. At the 
smallest value of the gas parameter (no^ = 10 “®), we 
find a good agreement between the VMG result and the 
GP prediction 


AEgp 


La;Ly-hcon 1 - 
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( 10 ) 


Significantly larger excitation energies are obtained at 
higher densities, showing that beyond mean-field ef¬ 
fects result in a more pronounced enhancement of the 
soliton energy Sill's] compared to the increase of the 
ground-state energy E\^u\ of the background homoge¬ 
neous gas [S^. From the slope of AE with respect to 
we extract the effective mass m* of the soliton shown in 
Fig. [5] as a function of na®. Despite the fact that the ex¬ 
citation energy becomes larger as the density increases, 
the effective mass per surface unit remains always close 
to the mean-field prediction = —4-|^ with a reduc- 
































tion of only ~10% at the largest density. We would like 
to point out that values of to* consistent with the ones 
reported in Fig. [5] are obtained from the linear depen¬ 
dence of the canonical momentum Pc in Eq. [7] on the 
velocity v of the soliton. 

In Fig. [5] we also show the mass mAN of missing par¬ 
ticles in the notch of a dark soliton per unit surface, ob¬ 
tained from the formula AN = {Ns — Njj)/L^Ly. We 
find that up to densities na^ ~ 10“^, the suppression of 
AN with respect to the GP value ANcp = is 

small and approximately equal to the reduction of the 
effective mass. At higher densities, the ratio AN/ ANqp 
drops down, following the filling of the soliton dip as 
shown in Fig.[TJ Due to the qualitatively similar behavior 
of TO* and mAN, their ratio remains surprisingly close to 
the mean-field prediction = 2 up to large values 

of the gas parameter {na^ ~ 10“^). Only at the largest 
density no? = 10“^ this ratio shoots up because of the 
large suppression in AN (see inset of Fig. [5]). It is worth 
stressing that the ratio m*/{mAN) is an important pa¬ 
rameter characterizing the dynamics of solitons as local¬ 
ized objects, being in particular related to the frequency 
of their oscillation in a harmonic confinement H- 

IV. CONCLUSIONS 

We have investigated the structural and energetics 
properties of soliton excitations in 3D using for the first 


time a fully microscopic, many-body approach and care¬ 
fully comparing our results with standard GP theory. We 
provide quantitative predictions for the effective mass 
and the mass of missing particles as a function of the 
gas parameter, showing that their ratio (that is the key 
parameter in the dynamics of the solitons in a harmonic 
confinement) is in good agreement with mean-field pre¬ 
diction even for regimes where the time-dependent Gross- 
Pitaevskii equation usually fails. 

In particular, the results at high density can be rele¬ 
vant for future experiments on the dynamics of solitons 
in strongly interacting composite bosons realized on the 
BEG side of a Feshbach resonance in a two component 
Fermi gas. At the moment, the main hindrance to the 
experimental determination of to*/(to AN) is the fast de¬ 
cay of the solitons, which limits the possibility of measur¬ 
ing the period of oscillation of these defects in harmoni¬ 
cally trapped gases [S^. By applying a stronger radial 
confinement to cigar-shaped configurations, it might be 
possible to increase the lifetime of solitons and thus to 
study their dynamics in regimes where the gas parameter 
is large. 
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